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1 Introduction 



By the exact discretization of an ordinary differential equation x = f{x), 
where x(t) G M^, we mean the difference equation Xn+i = F{Xn), where 
Xn G M^, such that X„ = a;(t„). Initial data coincide as well, because 
Xq = x{to). Exact discretizations have been first studied by Potts [28] and 
a detailed account of subsequent developments can be found in Agarwal's 
book [1]. It is worthwhile to point out that all linear ordinary differential 
equations with constant coefficients admit exact discretizations [H [28] . 

In this paper we discuss exact discretizations for the harmonic oscillator 
and its various extensions, including multidimensional, damped and inho- 
mogeneous cases. In other words, we consider the exact discretization of 
motions in a quadratic potential. 

Discretizations are closely related to numerical finite differnce schemes. 
Some numerical algoritms are known to integrate exactly the harmonic os- 
cillator equation, e.g., the so-called Gautschi-type methods [131 [ISl ED] and 
exponential integrators [221 ESI EZj- Recently, the problem of the exact dis- 
cretization of the harmonic oscillator equation with the constant force turns 
out to be very important for the construction of new numerical algorithms 
( "localy exact numerical schemes" [9l [HI [12] , see also Section 15. 4p and this 
is our main motivation for studying the subject of exact discretizations. 

2 Exact discretization of the one-dimensional 
harmonic oscillator 

It is well known [H [TUl [2B| that the discrete second-order linear equation 

+ iJ^Xn = (2.1) 



discretizes exactly the classical harmonic oscillator equation 

x + u;2x = 0, (2.2) 

where x = x{t) and the dot denotes the derivative with respect to t. In other 
words, any solution x„ of (12.11) can be expressed ), where the 

function t x{t) satisfies the equation (12.21) (the time step tn+i — tn = £ 



2 



is constant). Defining Vn '■= x{en) we derive the exact discretization of the 
velocity x(t) 

Xn+l — Xn COS US 

= . , 2.3 

UJ 

compare [TD] (see also below, where more general cases will be presented in 
detail). 

2.1 Equivalent forms of the exact discretization 

The exact discretization of the harmonic oscillator equation can be repre- 
sented in several equivalent forms which are useful in various contexts. The 
simplest form, following immediately from (12.11) . is 

Xn+l — 2 COS cue Xn + Xn-l = . (2.4) 

Equations similar to (12.41) can be found in [16] (Section 1.7, Examples 2 and 
5) without mentioning their relations with exact discretizations. 

Proposition 1. The equation ^2. ij) can he represented in the following equiv- 
alent form: 



Xn+l ~ COS CUE \ I Xn — Xn-l COS UJE 



smcje: 



coscje 

+UJ^Xn-l = . (2.5) 



Proof: A simple straightforward calculation. □ 



Corollary 2. Defining a new difference operator 
T — cos ojs 



^ ^ sin OJS 



(2.6) 



where T is the shift operator (i.e., {Tx)n = Xn+i), we can rewrite 5\) as 
A^x + cj^x = . (2.7) 
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Note that in the hmit e ^ the equation fl2.1l) assumes the form of the 
symmetric Euler finite difference scheme 



1 , , .2 



2 + Oo'Xn = , (2.8) 

while the equation (12. 5p apparently tends to the forward Euler finite differ- 
ence scheme: 

Xn+l Xn \ I Xfi Xn—\ 



e J \ e 



^ +u'xn-i = Q (2.9) 

(when cosue is approximated by 1). However, performing this limit with full 
care (i.e., leaving the second term \uj'^s'^ in the Taylor expansion of cosujs) 
we also get (12.81) . 

The harmonic oscillator equation (12.21) can be represented in the Hamil- 
tonian form 

mx = p, p = —kx , (2-10) 
where a;^ = — . 



Proposition 3. The exact discretization of Il2.1(j\) is given by 

m{Xn+l - Xn) Pn+l + Pn 



5 2 

Pn+l ~ Pn k(^Xn+l + X^) 



f2.111 



6 2 

where 

(5 = - tan — . (2.12) 
uj 2 

Proof: It is enough to show that the system ()2.1ip . where pn = rnvn, is equivalent 
to the system (j2.3p . (|2.4p . Solving ()2.1ip with respect to p„ and Pn+i> and then 
shifting index n + 1 ^ n in the first resulting equation, we get 

m k6 

Pn+l = -ri^n+l - Xn) -r{Xn+l + Xn) , 

4 

Pn = ^{Xn+l - Xn) + —{Xn+1+ Xn) , (2.13) 
4 

_ m \ ( \ 

Pn — ~r\Xn Xn—l) ~~r\Xn ~l~ Xn—l) ■ 
4 
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Eliminating p„ from the last two equations we get 

m k6\ f m k6\ / m k6\ , , 

- + - - 2 --- x„+ - + - = , (2.14) 



6 4 J \6 4 J " \6 4 

which is equivalent to (j2.4p if and only if 

which is satisfied by virtue of (j2.12p . Substituting (j2.12p and k = muj'^ into the 
middle equation of (I2.13P we get 

ma;(xn+i - Xncoswe) 

Pn = : , (2.16) 

smwe 



i.e., Pn = mvn, where Vn is given by (|2.3p . which ends the proof. □ 

Taking into account (I2.16P and the first equation of (12.131) we can solve 
the system (12. lip with respect to Xn+i,Pn+i obtaining 



Xn+i \ ^ f COS cue {mu) ^ sinus \ f \ (2 17) 

Pn+i J \ —ma; sin cue cos cue J \ Pr 

2.2 Discrete analogues of the energy integral 

The harmonic oscillator equation x + u'^x = has the following integral of 
motion 

^x^ + ^^^^a;^ = E = const , (2.18) 

where E can be interpreted as the energy per a mass unit. Trying to find 
a conservation law for its discrete analogue, (12. ip . one may consider the 
following quantities: 



eI^^ = x1j^_i — 2 cos ue XnXn+i + a; 



2 

n 1 



1 / ^n+l 

2 



^sin^ 



1 f Xn+l — XnCOSUJe\ 1 2 2 



(2.19) 



l.A ^2 1 

-(A^x„) + -uj~x: 
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Proposition 4. We assume that x„ satisfies \2.1\) . Then, E^^^^ = En for 
any k = 0,1,2 and any n G Z. Moreover, 




(2.20) 



Proof: The identities (j2.20p can be easily checked in a straightforward way. There- 
fore it is enough to show that En^ = E^}_-^. We compute: 

E^^ - E^l^ = {Xn+l - Xn-l){Xn+l + Xn-1 - 2 COS (J £ X„) , 

which vanishes provided that ()2.4p holds. □ 

Therefore all three quantities defined by formulas fl2.19p are integrals of 
motion, or, more precisely, they are different representations of the discrete 
analogue of the energy integral. In a sense these three forms of the energy 
integral correspond to the above three equivalent forms of the exact discrete 
harmonic oscillator equation. Namely, E^^^ can be obtained as a quadratic 
integral of motion for (12. 4p . Then, E^^^ is a natural guess in the spirit of 
Mickens' approach [24], strongly implied by the form (12. ip . Finally, E^"^^ is 
the best discrete analogue of the energy integral. Indeed, having the exact 
discretization we expect that there exists a discrete analogue of this conser- 
vation law, i.e.. 

En = \vl + \u^xl (2.21) 

where Vn has to be defined. Note that the invariant E'^'^'^ is of the form (12.211) 
provided that Vn is defined by (12. 3p . 

2.3 A family of geometric integrators 

In this subsection we consider a family of discrete systems containing the 
exact discretization of the harmonic oscillator. 

Proposition 5. The discrete map {xn,Pn) (yXn+i,Pn+i) , defined by: 

Xn+l - ^Xn + Xn-1 = , Pn = aX^+l - (3Xn , (2.22) 
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is both symplectic and energy-preserving for any choice ofa,(3,'y. If, more- 
over, the coefficients depend on the time step e, satisfying constraints 



a{-e) = -a{e) , 7(^) = ^^'^ .^j , (2.23) 
then the map is also symmetric (time-reversible). 

Proof: We compute Pn+i = a(7a^n+i — — Pxn+i- Hence, 

dXn+l A dpn+i = -adxn+1 A dXn , 
dxn A dpn = adxn A dXn+1 ■ 

Thus dxn+1 A dpn+1 = dxn A dpn which means that the map is symplectic (and 
volume-preserving) [15j. An integral of the discrete evolution can be easily found 
by multiplying the first equation of (j2.22p by Xn+i — Xn-i- Indeed, 

= {Xn+l - JXn + X„_i)(Xn+l - Xn-l) = X^+i - JXn+lXn + 7X„X„_i - xl_i , 

which means that 

x^+i - 7x„+ix„ + x^ = x^ - 7x„x„_i + xl_i = const . (2.24) 
In order to find conditions for time-reversibility [15], we rewrite ()2.22p as follows: 

H 1 

a a 

7/3-^ 7-^ 
where a = a(e), /? = /3(e) and 7 = 7(e). Computing we get: 



/ Xn+l 1 


1 = A{e) 1 


( Xn \ 


\ Pn+1 1 




[ Pn 1 



A{e) 



I Xn 1 


1 =^-^(^)| 


( Xn+l \ 






{ Pn+1 I 



A-\e) = 



' a a 

^-7/3 ! 



Time-reversibility means that A{—e) = A ^(e), i.e., 
a{—e) = —a{e) , 

7(-) - = , (2.25) 
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The first two equations yield (j2.23p . The third equation is identically satisfied 
provided that the first two equations hold. □ 



The maps defined by fl2.22p contain a large family of geometric numerical 
integrators for the classical harmonic equation (12.21) . Indeed, if 

2 — 'y(e) k 

\imea(e) = m , \ime3(e) = m , lim = — = uj'^ , (2.26) 

s^O e^O e^O m 

then the continuum limit of (I2.22p (divided by e^) is given by x + uo'^x = 
and p = mx. In particular, the exact discretization of the harmonic oscillator 
is characterized by 

TfiUJ 

aie) = , l3(e) = mujcot(uje) , 7(e) = 2cosco'£: . (2.27) 

sm Lue 

Corollary 6. The exact discrtization of the harmonic oscillator equation sat- 
isfies all assumptions of Proposition and, therefore, is symplectic, volume- 
preserving, energy-preserving and time-reversible. 



3 Harmonic oscillator with a constant driv- 
ing force 

In this section we consider the exact discretization of the classical harmonic 
oscillator equation with a constant force: 

X + cu'^x = g , (3.1) 

where u and g are constant. This problem deserves special attention in the 
context of some numerical applications, see Section [51 

3.1 Variable time step 

The general exact solution of the equation (13. ip and its derivative read 

x{t) = A cos cut + -B sin cut + ^ , 

^ (3.2) 



x{t) = ooB cos cut — 00 A sin cut . 
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We consider the exact discretization 

(3.3) 

Vn '■= (^B COS utn — ujA sin utn ■ 



Proposition 7. The exact discretization (with a variable time step) of the 
equation ( 15'. 1\) is given by 

O'n-lXn+l ~ bjiXn + 0-nXn-l = — n {O'n + C^n-l ~ ^n) ; (3-4) 

where an = sinu;£:„, 6„ = sinco'(£:„ + En-i) and En is a prescribed (variable) 
time step. 



Proof: Assuming the variable time step, tn+i — tn = £n, and using ()3.3p . we have 

(3.5) 



Xn+l = A COsiiOtn + UJEn) + B Smiujtn + OJEn) + ^ 



Xn-l = Acosiijjtn - OJEn-l) + B Smiujin - OJEn-l) + ^ • 

Therefore 

Xn+l = {Xn %) COS UJEn + sin (S COS wt^ - Asiniotn) + 

(3.6) 

Xn-i = (xn ^)coscje„__i - sin u;e„_i (S COS a;t„ - Asina;t„) H ^, 

and adding the first equation muhiphed by sina;e„_i to the second equation mul- 
tiphed by sina;e„ we get 

Xn+l - A) sina;en-i+ (xn-i - A) sina;e„ = {xn- -^)sin{uj£n + u;£n-i) , 

which is equivalent to ()3.4p . □ 

Comparing (13. 3p and the first equation of (13.61) . we obtain 

_ X„+i - Xn COS UJEn - ^ (1 - COS U)En) 

which is the exact discretization of the velocity x. 



Proposition 8. The exact discretization of the harmonic oscillator with a 
constant force 1\) is given by 

Xn+l \ f COSUJSn OJ'"^ SUmJEn \ f \ 9 ( I — COS US r. 



Vn+1 J \ —ujsm.ujen COS uJEn J \ '^n I oj"^ \ a;smci;£: 



(3.8) 



where Vn is the exact discretization of the velocity x and the time step Sn is 
variable. 



Proof: From (|3.3p we derive: 
Xn+l = {A COS UJtn + B sin Ujtn) COS UJEn + {B COS LUtn — A sin LOtn) sin LOSn + oj~'^g, 

Vn+i = —{Aijo cos ujtn + Buj s'm ujtn) sin we„ + {Buj cos ujtn — Auj sinwt^) cos wen, 
Hence 

Xn^\ — Lu^'^g = {xn — i^^^^) COS a;e„ + Vn^~^ sma;e„ , 

(3.9) 

Vn+l = -{Xn-UJ ^g)uJ Sill UJEn + Vn COS UJEn , 

which is equivalent to (|3.8p . □ 

Proposition 9. The exact discretization of the equation ^3. 1\) can be repre- 
sented in the following form, equivalent to \3.3^) : 



T = o '^^'^+1 + ' 

On ^ 

Vn+1 — Vn 1 2 / , \ , 

- --UJ [Xn+l + Xn)+ g 



(3.10) 



5n 2 

where 

Sn = -tan^^ . (3.11) 
Proof: Analogical to the proof of Proposition [3l □ 
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3.2 Constant time-step 

The equation (13.41) simplifies in the case e„ = e = const. It is convenient to 
use the identity 

a„ + a„_i - 6„ = 4 sm ^ sm — sm — ^ — . (3.12) 



Corollary 10. The equation 

Xn+i-'2xnCosue + Xn-i = '^sm^^ (3.13) 
is the exact discretization of \3.1\) with a constant time step e. 



Proposition 11. The equation Ii3.13\) can be rewritten in the following equiv- 
alent forms: 



+ uj'xn = g, (3.14) 



: \- 00 Xn-1 = ■, (3.15) 

smcjg " ^ COS — 

UJ 2 



ff„ — tan ^) — (t;„_i — tan ^) coscje „ 

\ZJ ^ + = 9 , (3.16) 

UJ 

where f„ is defined by /{2.3\) . 

3.3 Discrete analogues of the energy integral 

The energy integral for the equation (13. ip is given by 

E =^x'^ + ^u^x^ - gx . (3.17) 
Proposition 12. Assuming that Xn satisfies Ii3.13\) we define: 

En'^ = - 2C0SU;£ X„X„+i + X^ - {Xn+Xn+l) 9, (3.18) 
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■■= 2 [-J—If ) +2^^n^"+i 2 9, (3.19) 

1 / - x.cos^e \ ^ 1^ _ x^^ + i^ 

2 V / 2 " 2 cos2 ^ ^ ' ^ ^ 

n 2 V ^li^ u; 2/2 " y " ^ ^ 

Then, -E^^i = -E'i'^^ /or any A; = 0, 1, 2, 3 and any n G Z. Moreover, 

2 V2smf " ' 
^...(l,,.,,^)4..4(_^)^^.o, (3.22) 

= Sf +(5ta„if)' , 

Proof: The identities p.22p can be easily verified. Therefore it is enough to show 
that E^^ = eI^I^. We compute: 

- E'lf^i = (Xn+l - Xn~l) (^Xn+l + Xn-l - 2 COS - ^ siu^ , 

which vanishes provided that (j3.13p holds. □ 



Therefore all four quantities defined in Proposition [T^ can be interpreted 
as discrete analogues of the energy integral. They correspond to different 
equivalent forms of the exact discrete harmonic oscillator equation with a 
constant driving force. In particular, E^^^ seems to be natural in the frame- 
work of Mickens' approach [23]. However, the best discrete analogue of the 
energy integral is given by E^^\ Indeed, E^^^ is of the form (13.171) evaluated 
at X = Xn and x = Vn, where 



x„+i - x„ cos cue g ue 
Vn = ■■ tan — , [6.2.6) 

" sma;£ 2 



compare (13. 7p . 
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3.4 Damped harmonic oscillator 

The damped harmonic oscillator with a constant driving force 

X = —UqX — 2jx — g (3.24) 
can be reduced to the harmonic oscillator without damping. 
Proposition 13. The transformation 

X = e^* {x - Xe) , (3.25) 
where Xe = —u^'^g, reduces iS.24\) to the harmonic oscillator equation 



X + cj^X = 0, uj^ = uj^Q-Y- (3.26) 

Proof: Straightforward computation. □ 

Denoting p = mx and P = mX we express P in terms of x, p, namely: 

P = e^^ {p + m-f{x - Xe)) . (3.27) 

Therefore, the exact discretization of the equation fl3.24p in terms of variables 
X, P is given by 



Xn+i \ _ f cosue (mu) -"^ sin cue \ / X„ 
Pn+i J \ —muj sin ue cos ue y y P„ 

(compare fl2.17p ). Finally, substituting 

Xn e^ " i^Xfi Xe) , 

Pn = e"^*" {p„ + m-f{Xn - Xe)) 

into fl3.28p we obtain the final corollary. 



(3.28) 



(3.29) 



Corollary 14. The exact discretization of the damped harmonic oscillator 
equation with the constant driving force ^3.24^ is given by 

Xn+l = Xe + e~^'^ {{Xn - Xe) {cOSUS + ^ siu + p„) , 

Pn+i = e '^^ [pn (cos we - ^ smue) - ( + ^ j (x„ - Xe) smue 
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4 Exact discretization of the multidimensional 
harmonic oscillator equation 

In this section we consider the equation 

^+n'x = m, (4.1) 

where x = x(t) G M", is a given invertible constant n x n matrix and / = 
/(t) G M" is a given driving force. This is a natural n-dimensional extension 
of the harmonic oscillator equation. Indeed, if / = and Q"^ has n pairwise 
different eigenvalues, then the equation (14.20 describes n independent one- 
dimensional harmonic oscillators. 

Here, in contrast to previous sections, we try to use the variable time-step 
whenever possible. 

4.1 Multidimensional harmonic oscillator 

We begin with the case f(t) = 0, i.e., 

§ + tf.^O. (4.2) 

It is convenient to represent (14.21) as the following first order system 

X = V , V = — . (4.3) 
The general exact solution is given by 

x{t) = e*™ci + e-*^*C2 , 
v{t) = in (e^™ci - e-^^ca) , 



(4.4) 



where Ci, C2 are constant n- vectors. 

The exact discretization of the equation (14. 2 p is given hy Xn '■= x{tn) 
where we use ( 14.41) . Therefore 



(4.5) 
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Hence 

e**"% = ^ K - zl^" V) , e-**"% = ^ (x„ + zr]- V) , (4.6) 
and, denoting En = tn+i — tn (the time step, in general variable), we have 
Xn+i — a Ci -\- e C2 , 

Substituting (14 .Gp into (14.71) we obtain: 



(4.7) 



Xn+l \ _ ( COS QSn ^ siu \ ( 

Vn+i J I — f2sinfi£:„ cosfie^ ) \ v„ 



(4.8) 



where we use a natural matrix notation (the matrix entries and vector com- 
ponents are n x n matrices). In particular, we have: 

Vn = QisinQen)'^ {xn+i - (cosil£:„)x„) . (4.9) 



The matrix on the right-hand side of (14.81) is an even function of Q, i.e., this 
matrix can be expressed as a series in terms of Q"^ (without explicit knowledge 

of n). 

Corollary 15. The exact discretization of the system lji4-3\ ) is given by (^7^, 
where En is an arbitrary variable time step (in particular, we can take En = 
E = const ). 

Below we present three propositions. The proofs are omitted because in 
Section W72\ we will prove more general theorems. 

Proposition 16. The exact discretization ( [^.<j^ of the multi- dimensional 
harmonic oscillator equation Iji4-S\ ) is equivalent to the system 

{Xn+l -Xn) = I {Vn+l + Vn) , 

(4.10) 

5„ ^ {Vn+l - Vn) = {Xn+1 + Xn) 



where 



2fi-itan^. (4.11) 



15 



Proposition 17. // the time step is constant = ^j, then the exact dis- 
cretization ( [^.<S] j of the multi- dimensional harmonic oscillator equation Iji4-^ 
is equivalent to the system 

Xn+l - 2(cOsfi£:) Xn + Xn-1 = , 

(4.12) 

Vn = ^{smVle) (x„+i — (cosf2£:)a;„) . 
Proposition 18. // f„ is defined by Iji4-^ o.iT'd = ^, then 

In:=^\Vnf + ^{Xn\^^Xn) (4.13) 

is an integral of motion of the discrete multidimensional harmonic oscillator 
equations ^.10^ (i.e., = In). 

Here (and below) (01 ip) denotes a scalar product in M", |0p = {<f)\ <f)) 
and is the transpose of fl. We recall that 

(M0 1 ^) = (0 1 M^^p) and {(j)\^) = {^p\ 0) (4.14) 

for any matrix M and any vectors 0, "^/^ G M". Moreover, it is worthwhile to 
point out the obvious fact that Q, 6n, sinQsn (and other analytic functions 
of Q) pairwise commute. We use frequently this property. 

4.2 Multidimensional harmonic oscillator with a con- 
stant driving force 

We consider the equation 

^ + n^x = a, (4.15) 

where x = x{t) G and f2 is a given invertible constant n x n matrix and 
a = const e M". It is convenient to represent (I4.15P as the following first 
order system 

X = V , V = —Q'^x + a . (4.16) 
Its general exact solution is given by 

x{t) = e*™ci + e-*™C2 + n-^a , 

(4.17) 

v{t) = in (e*^*ci - e-*^*C2) , 

16 



where Ci, C2 are constant n- vectors. 

The exact discretization of the equation fl4.15l) is given by := 
where we use f l4.17p . Therefore 

(4.18) 

Vn = (e**"^Ci — e **"^C2) . 

Hence 

(4.19) 

1 

Denoting = — t„ (the time step, in general variable) and evaluating 
f l4.18p at n + 1, we obtain 

Xn+i = e**"^+*^"^ci + e"**"™"*^"^c2 + r^^^a , 

(4.20) 



Proposition 19. The exact discretization of the system Iji4-i6\ ) is given by 



Xn+l \ _ f COS ile-n f2 ^ siu \ f Xn \ , f 2i7 ^ siu^ ^ 



a 



Vn+i J V sin flSn COS f2£„ J \ "^"^ J \ ^ sin Qsn a 

(4.21) 

where En is an arbitrary variable time step (in particular, we can take En = 
E = const ). 

Proof: It is enough to substitute (|il9]) into ([OO]) . □ 
In particular, 

Vn = n(sin ri&„)~"^ {xn+i — (cosri£„)x„) — tan — — ^ a . (4.22) 
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Proposition 20. The exact discretization ^ ■ 21\ ) of the multi- dimensional 
harmonic oscillator equation with the constant driving force ( [^.16] ] is equiv- 
alent to the system 



5„ ^ {Xn+l - Xn) = \ {Vn+l + Vn) , 

{Vn+1 - Vn) = -l^^ {Xn+1 + Xn) + a 

where 6n is given by Ili4-ii\ )- 



(4.23) 



Proof: From (j4.2ip we compute: 
Xn+i + Xn = {I + COS Q,en)xn + {^~^ siii Oe„)t;„ + (1 — COS ile„)0~^a , 

Xn+i — Xn = —{I — cosQ,en)xn + (O^"*^ sin + (1 — cosf]e„)^~^a 

Vn+l + Vn = —{^ sin Qen)Xn + (1 + COS Qen)vn + siii QSn a , 

Vn+1 — Vn = — (il sin r2e„)xn — (1 — COS ^len)vn + 0~-'^sinOe„ a . 
Hence 

. ^£n f I \ I I \ no— 1 ■ 

SZsm-^ (X„+1 + Xn) + COS (t^n+l - Vn) = 2iZ ^m a , 



(4.24) 



(4.25) 



Vl cos ^ (xn+i - 2;„) = sin ^ (i;„+i + ?;„) , 
which is equivalent to (j4.23p . □ 

Proposition 21. // the time step is constant (sn = then the exact dis- 
cretization fi4-21^ of the multidimensional harmonic oscillator with a constant 
force can be rewritten in the following equivalent form 



(4.26) 



2Q,~^ sin — j a , 

Vn = ^l{sm^le)~^ (xn+i — (cosf2e)x„) — r^^^tan^ a . 
Another convenient expression for the velocity f„ is given by: 

Vn = ^fi(sinfie)~^ (x„+i - x„„i) . (4.27) 
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Proof: From (j4.2ip we get: 

^}^^{sinUe)vn = {xn+i — {cosQe)xn) — il^^tan^ a , 

Vn = sin ile + cos Q,e v^-i + sin J7e a , (4.28) 

Vn-i = r2(sin r^e)"-^ (x„ — (cos ile)a;„_i) — tan ^ a . 

Thus we derived the formula (j4.26p for Ehminating Vn and Vn-i from the 
system (|4.28p we get the first equation of ()4.26p . In order to obtain the formula 
(14271) we evaluate (l4T8]) at n ± 1 

= e**"^±'^"^ci + e-^*"f^*T*^"f^c2 + Vl-'^a , (4.29) 

and compute 

^ (x„+i - x„_i) = ie*"^(sine„Jl)ci - ie-^*"^(sine„0)c2 . (4.30) 
Taking into account the second equation of ()4.18p we complete the proof. □ 



Proposition 22. // Q = Q and f„ is defined by 11^4 -2 2^ , th 



len 



In ■■= \ + ^(a^nl ^"^Xn) - (o | X„) (4.31) 

is an integral of motion (i.e., In+i = In) of the discrete multidimensional 
harmonic oscillator equations \4-2'^- 



Proof: From (j4.23p we have 

Vn+l +Vn = 26n^ (x„+i - X„) , 

1 2 (^-^^^ 

Vn+l -Vn = -2^n^ {Xn+1 + Xn) + . 

In order to prove In+i = In it is enough to multiply the above equations side 
by side (using the scalar product and its properties). = O obviously implies 
= ^n- Taking also into account ()4.14p . we verify that 

I 6n^'^Xn) = {Sn^Xn \ 6n^'^Xn+l) ■ (4.33) 

Therefore, multiplying ()4.32p side by side, we get 

= I dn^'^Xn+l)-{5~^Xn \ 5n^'^Xn)+2{Xn+l - X„ | o) . (4.34) 

Hence, applying once more (|4.14p . we get In+i = In- ^ 
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4.3 Multidimensional harmonic oscillator with a poly- 
nomial driving force 

In this section we consider the harmonic oscillator perturbed by an arbitrary 
polynomial driving force: 

N 

x + n^x = f{t), f{t) = Y,Ckt\ (4.35) 

fc=0 

where x = x{t) G M", f{t) G M", ^ and is a constant invertible 
matrix. 

The key observation is that the general solution to fl4.35p can be repre- 
sented in the form 

x{t) = (cos nt)A + (sin nt)B + $(t) , (4.36) 

where where A, S G M" and 

oo 

<l>{t) ■= ^(-f]-2)^'/(2fc)(i) . (4.37) 

k=0 

One can easily see that the infinite sum contains only a finite number, [A^/2], 
of non-zero terms and $(t) is a polynomial of A^th order. Moreover, 

v{t) = x{t) = fi(cos nt)B - Q{sm Qt)A + $(t) . (4.38) 

Proposition 23. The exact discretization of lji4-35\ ) is given by 

Xn+l - 2 COS Qe Xn + Xn-l = $„+l - 2(cOS Qe)^n + , 

(4.39) 

Vn = ^Q{smQe)'^ (x„+i - - + <l>„_i) + $„ , 
where = $(ne), = ^{ne) and $(t) is defined by 

Proof: Computing x{t it e) and expressing A, B in terms of x, v, <I> we get: 
x{t±e) = {cosne){{x{t) - ^{t)) ± n-\smne){v{t) - ^{t)) + ^{x ± e), 
x{t + e) + x{t -e)= 2(cos Oe)(x(t) - ^{t)) + $(t + e) + $(t - e) , (4.40) 
x(t + e) - x{t -£)= 21^-1 (sin 17e)(i;(t) - 6(t)) + ^{t + e) - ^{t - e) . 
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Now, identifying t — > ne, x{ne) x„, $(ne) we obtain the exact discretiza- 

tion of the perturbed harmonic oscillator equation (j4.39p . □ 

In the particular case f{t) = a = const, we have 
$(t) = Q-y = Q-^a , 

(4.41) 

- 2cosQe^ri + ^n-i = 2(1 - cosQe)Q . 
Therefore, in this case (14.391) reduces to (I4.26p . 

Analogical considerations can be made for any other perturbation (driving 
force) f = f{t) such that the series 

oo 
fc=0 

is finite or summable. For instance, in some cases we get geometric series. If 
/(t) = e°*/o, (a = const G M), then 

oo 

$(t) = ^(-a2fi-2)'=e"Vo = {^^ + e"Vo , (4.43) 

where I is the unit matrix in M". If f{t) = /osintut, {uj = const G M), then 

oo 

$(t) = ^]-2^(cu2^]-2)'^(sincut)/o = {Q^ -ujHy'e'^'fo • (4.44) 

k=0 

Using the linearity of (14.371) one can extend these results on finite Fourier 
series, finite linear combination of exponential functions or some specific com- 
binations of polynomials, exponentials and triginometric functions. Similar 
results were first obtained (in a different way) a long time ago [21 

5 Connections with numerical methods 

In order to produce an exact discretization we have to know corresponding 
exact solutions and in such cases numerical schemes seem to be of a little 
practical use. It turns out, however, that exact discretizations can be applied 
to construct good numerical schemes which are exact in some particular cases. 
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5.1 Gautschi-type methods 

The exact discretization of the harmonic oscillator is of great importance 
for a large class of numerical methods (Gautschi-type methods) |13l [18]. In 
order to illustrate this idea we consider the following class of equations: 



where a; G M", a; G M and the nonlinear term g{x) is assumed to be much 
smaller than the linear term uj'^x. The simplest Gautschi method is defined 
as 



One can easily see that for g{x) = const the Gautschi method [13] reduces 
to the exact discretization of the harmonic oscillator with a constant driving 
force, compare ( 13.13^ . In other cases the Gautschi method ceased to be exact 
but is still considered as one of the best for oscillatory problems with constant 
high frequencies [I5l [18] . 

5.2 Exponential integrators 

Exponential integrators [25l [27] are also concerned with exact discretizations. 
We confine ourselves to equations of the form 



where 7/ G M"-, L is a constant matrix and g{y) is a small nonlinear term. We 
consider three exponential integrators which are extensions of the classical 
Euler method: the explicit Euler-Lawson scheme 



X = —uP'x + g{x^ 



(5.1) 




(5.2) 



y = Ly + g{y) , 



(5.3) 



Vn+i = e'^ {jjn + eg{yn)) 



(5.4) 



the implicit Euler-Lawson scheme 



Vn+i = e^^yn + eg{yn+i) , 



(5.5) 



and the exponential Euler method: 



yn+i = e'^yn + eipi{sL)g{yn) , 



(5.6) 
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where ipi{z) := (e^ — l)/z. The equation fl3.1l) obviously can be represented 
in the form fl5.3p : 



±f x\ _ f 1\ f x\ f 
dt[v -[-uj' J U j + U 



(5.7) 



where g = const. Applying the above exponential integrators to (15. 7p we get, 
respectively, 

Xn+i \ — f cosue u'^sinue \ f \ _^ eg f sin cue \ 
Vn+i J \ — djsincue cos lue ) \ '^n ) V ^ COS cue: J ' 

Xn+i cosujs LU'^sinuje W a;„ A _^ / \ 

Vn+i ) V sin cue co&uje ) \ Vn ) \ eg J ' 

Xn+i — uj"'^g \ ( COS cue uj~^ smuje \ ( Xn — uj^'^g 



Vn+i J \ —UJ sm one cos ue 

We see that in the case g = (the harmonic oscillator without the driving 
force) all three integrators yield the exact discretization ( 12.17p . In the case 
g ^ only the exponential Euler scheme (15.61) is exact, compare (13. 8p . 

5.3 Application to the Kepler motion 

Sometimes considered equations reduce to the harmonic oscillator (or contain 
the harmonic oscillator equation). Then we may take advantage of using the 
exact discretization of the harmonic oscillator. Here we shortly describe the 
case of the classical Kepler problem 

mr = J , r G , r = |r| , (5.9) 

where m,k are constant. It is well known (see, for instance, [I5l|29|), that 
exact trajectories are given by r = r (</?), where ip is the angle swept out by 
r and 1/r solves the so called Binet equation (equivalent to the harmonic 
oscillator equation with the constant driving force): 

(Pu km 1 

where L = |r xmr| = const (the angular momentum) is an integral of motion. 
Using the exact discretization of the harmonic oscillator equation f l5.10p we 
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get a numerical scheme such that (pn+i — = const (but the time step is 
variable) [3. 

It is also well known [20j that the Kepler problem can be reduced to the 
4-dimensional harmonic oscillator equation 

^-|g = o, g = Q(s)GR\ (5.11) 

where E is the energy integral of the considered Kepler motion. Having Q 
and s we can find r (explicitly expressed in terms of Q) and t (integrating 
the equation i! = IQp), for details and precise formulas see, for instance, 

[laisniiss]. 

However, the possibility of using the exact discretization has been often 
overlooked. First, Minesaki and Nakamura [26j discretized the harmonic 
oscillator equation (15.111) by the discrete gradient method [HI [21]. The 
obtained discretization was quite good (the exact trajectories were preserved) 
but not exact. Then, Kozlov [19] succeeded to get the exact discretization 
by summing up some infinite series. Earlier papers [31 [30] , where the exact 
discretization of the harmonic oscillator was taken into account, now seem to 
be forgotten (as far as geometric numerical integration is concerned). More 
detailed discussion of the exact discretization of the Kepler motion based on 
the Kustaanheimo-Stiefel theory can be found in [8]. 



5.4 Locally exact modification of the discrete gradient 
scheme 

Locally exact numerical schemes generalize the well known concept of the 
small oscillations approximation [9l [TH [12] . If a numerical method contains 
some parameters, we modify these parameters so that the resulting scheme 
is exact for small oscillations [11] or, in more general setting, we require that 
the resulting scheme is exact for the linearization of the considered system 
[9], [12] . As an illustrative example we consider the system 

i] = -$'(x) , V = X , (5.12) 

where $ = $(a;) G M is a given function. The case $(x) = ^uj'^x^ + gx 
corresponds to the harmonic oscillator with a constant driving force. The 
energy conservation law reads 

]-v^ + ^{x) = const . (5.13) 
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The discrete gradient method [TH [2T1 [23] for fl5.12p yields: 

Vn+l - Vn _ ^jXn+l) - $(x„) 

(5.14) 

- {Vn+1 + Vn) = . 

^ s 
We consider the following extension of the discrete gradient scheme: 

Vn+l - Vn _ - ^{Xn) 

(5.15) 

- [Vn+l + Vn) = T , 

2 dn 

where Sn is an arbitrary positive function of en,Xn,Vn,Xn+i,Vn+i etc. The 
system f l5.15p is a consistent approximation of (15.121) if we add the natural 
consistency condition 

lim— = 1. (5.16) 

e^O e 



Proposition 24. The numerical scheme Ii5.15\) preserves exactly the energy 
integral (for any positive function 6n)- 

Proof: We multiply side by side both equations of (jS.lSp obtaining: 

^vl^^ + $(xn+i) = + $(x„) , (5.17) 

which ends the proof. □ 

In order to find the most accurate discretizations among (15.151) we lin- 
earize (15.151) around x = x (and x will be specified below): 

Cn+l 1 / . N 

7 = o iPn+1 + Pn) ■ 

On 2 

(5.18) 

= -$'(x) - ^$"(x) (e„ + M , 

On 2 
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where ^„ := Xn — x and = Xn+i — x. The resulting system (15.181) is 
equivalent to the harmonic oscillator equation with a constant driving force, 
compare fl3.10p . provided that we identify p„ = f„ and require 

g = , a;2 = V"{x) , (5„ = - tan — (5.19) 

U! 2 

(note that the consistency condition (15.161) is obviously satisfied). Choosing 
X = Xn or X = |(x„ + Xn+i) (thus changing x at every step) we obtain two 
very accurate ( "locally exact" ) numerical schemes, compare [12] . 



5.5 The wave equation 

In this section we are going to show that the exact discretization of the har- 
monic oscillator equation can be applied also in the case of partial differential 
equations. We consider the linearized wave equation 

utt = u,^^ -a^u . (5.20) 

Folowing Bridges and Reich ([3], Section 4.2) we perform the Fourier trans- 
formation, substituting u{x,t) = u{t)e^^^ . Hence 

= -uj\ , (5.21) 

where a;^ = k"^ + a?. Instead of using standard methods (like implicit mid- 
point/trapezoidal schemes, see [5]) we take advantage of using the exact 
discretization: 

M"+^-2(coscuAt)M" + u"-i = . (5.22) 

Thus the "numerical frequency" of the solution (see [5]) coincides with 
uj. It leads to many important improvements. For instance, the group veloc- 
ity becomes exact and our numerical scheme reproduces exactly the energy 
transport in wave packets. 



6 Conclusions and future directions 

The exact discretization of linear ordinary differential equations is a subject 
rather well known but a little bit forgotten. In this paper we have shown a 
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large number of useful properties and equivalent formulations of the exact 
discretizations of harmonic oscillator equation and its variuos extensions and 
generalizations. Among these results it is worthwhile to distinguish Propo- 
sitions [16] and [201 with formulas fl4.10p and fl4.23p . and Proposition [22] (the 
energy conservation law in a quite general case). 

Our main motivation for careful studies of exact discretizations consists 
in potential numerical applications. This point is explained in Section [3 In 
particular, a new concept for modifying numerical schemes is presented in 
Section 15.41 We point out that this derivation of the locally exact discrete 
gradient method is short and much simpler than the presentation given in 
our paper [12]. This simplification is due to the direct use of the formula 
f l3.10p which is especially suitable for any applications related to the discrete 
gradient method. Other results of this paper turn out to be very helpful in 
deriving and studying other locally exact numerical schemes, to be published 
soon [9] 

A promising future direction is associated with the time scale approach 
[U, [T7|. Its main goal is to unify differential and difference calculus. Exact 
discretizations also connect smooth and discrete case but in an apparently 
different way. We plan to find links between these two approaches. A nat- 
ural possibility is to generalize the notion of the delta derivative using new 
difference operators, e.g., the difference operator A^, see (12. 6p . Exact dis- 
cretizations of the harmonic oscillator equation can be also applied in the case 
of operators on Banach spaces (see [B]) and for some PDEs (see Section [S]). 
In the near future we plan to study in more detail the exact discretization of 
the wave equation. 
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